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' We extend the study of the light hadron spectrum and the quark mass in two-flavor QCD to 

smaller sea quark mass, corresponding to mps/mv = 0.60-0.35. Numerical simulations are carried 
out using the RG-improved gauge action and the meanfield-improved clover quark action at /3 = 1.8 
£Nj ' (a = 0.2 fm from p meson mass). We observe that the light hadron spectrum for small sea quark 

mass does not follow the expectation from chiral extrapolations with quadratic functions made 
from the region of mps/mv = 0.80-0.55. Whereas fits with either polynomial or continuum chiral 
jyrj ' perturbation theory (ChPT) fails, the Wilson ChPT (WChPT) that includes a 2 effects associated 

with explicit chiral symmetry breaking successfully fits the whole data: In particular, WChPT 
correctly predicts the light quark mass spectrum from simulations for medium heavy quark mass, 
such as mps/mv > 0.5. Reanalyzing the previous data with the use of WChPT, we find the mean 
up and down quark mass being smaller than the previous result from quadratic chiral extrapolation 
by approximately 10%, mjff = 2 GeV) = 3.11(17) [MeV] in the continuum limit. 
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Recent years have witnessed steady progress in the lattice QCD calculation of the light hadron spectrum [lj. 
In the quenched approximation ignoring quark vacuum polarization effects, well-controlled chiral and continuum 
extrapolations enabled a calculation of hadron masses with an accuracy of 0.5-3% 0. At the same time the study 
established a systematic deviation of the quenched light hadron spectrum from experiment by approximately 10%. 
We then have made an attempt of full QCD calculation that allows chiral and continuum extrapolations within a 
• consistent set of simulations @ • The deviations from experiment in the light hadron spectrum are significantly reduced 
and the light quark mass decreases by about 25% with the inclusion of dynamical u and d quarks. With currently 
available computer power and simulation algorithms, however, the sea quark mass that can be explored is far from 
the physical value and a long chiral extrapolation is involved to get to the physical u and d quark mass. 

An attempt has been made to push down the simulation to a small quark mass corresponding to mps/my rs 
0.3 in full QCD with the Kogut-Susskind(staggered)-type quark action |4j|. The staggered action, however, poses 
a problem of flavor mixing, which would modify the hadron spectrum and its quark mass dependence near the 
chiral limit. The staggered action also suffers from ambiguities in hadron operators and has a potential problem 
of non-locality. The Wilson-type quark actions have the advantage of simplicity: they are local and respect flavor 
symmetry, but a larger computational cost limits the simulations to relatively large quark masses corresponding to 
mps/mv > 0.6 0,01111 00 El El El- An important problem is to examine whether chiral extrapolations from such 
a quark mass range lead to results viable in the chiral limit. 

Chiral extrapolations are usually made with polynomials in the quark mass. The problem is that they are not 
consistent with the logarithmic singularity expected in the chiral limit. In reality, the physical quarks are not exactly 
massless and hence the polynomial extrapolation should in principle work. However, increasingly higher orders are 
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needed should one wish to increase the accuracy of the extrapolation. It is compelling to estimate the systematic 
errors due to higher order contributions when the data are extrapolated using a low-order polynomial. 

An alternative choice for chiral extrapolations is to incorporate chiral perturbation theory (ChPT) 0] . The present 
lattice data, however, are not quite consistent with the ChPT predictions. The high-statistics JLQCD simulation 
of two-flavor full QCD, using the plaquette gauge action and the 0(a)-improved Wilson quark action at j3 = 5.2 
(a = 0.0887(11) fin; the spatial size L ~ 1.06-1.77 fm), shows no signature for the logarithmic singularity in the 
pion mass and pion decay constant 11]. A possible reason for the failure to find the chiral logarithm is that sea 
quark masses, corresponding to mps/my = 0.8-0.6, are too large. Higher order corrections of ChPT may have to be 
included to describe the data, as suggested from a partially quenched analysis, which shows that mps/my = 0.4-0.3 
is required for the convergence of one- loop formula |13L Il4| . Another possibility is explicit chiral symmetry breaking 
of the Wilson quark actions that may invalidate the ChPT formulae. Modifications due to finite lattice spacings may 
be needed for an analysis of data obtained on a coarse lattice. 

Recently studies were made to adapt ChPT to the Wilson-type fermion at finite lattice spacings (WChPT) 0, 
HH H21 HI] > with subtle differences in the order counting, and hence the resulting formulae for observables, among the 
authors. The work assumes the O(a) chiral symmetry breaking effects being smaller than those from the quark 
mass, and only the effects linear in lattice spacing are retained in the chiral Lagrangian. This contrasts to the authors 
of Refs. [TtI lisj who include the 0(a 2 ) effects in the chiral Lagrangian, however, with different order countings. In 
Ref. p"?! the O(a) terms are treated as being comparable to the quark mass term while the 0(a 2 ) terms are assumed 
to be subleading: in this case, 0(a) effects are essentially absorbed into the redefinition of the quark mass in the 
one-loop formulae and the 0(a 2 ) terms provide additional counter terms. In Ref. 18], on the other hand, the terms of 
0(a 2 ) are kept at the leading order, because the existence of parity-broken phase and vanishing of pion mass depend 
on them in a critical way |15| . The coefficients of chiral logarithm terms receive 0(a) contributions, and hence the 
logarithmic chiral behavior is modified at a finite lattice spacing. Similar attempts to include the 0(a 2 ) flavor mixing 
for the staggered-type quark action were made in Refs. (19t l20l l2lj . 

The qq+q collaboration j^] applied the one- loop ChPT and WChPT with the prescription of Refs. [l^.ll7l| to their 
data obtained at mps/my = 0.9-0.5. Their simulations were made at coarse lattices of a = 0.19 fm (f3 — 5.1) and 
0.28 fm (f3 = 4.68) using the plaquette gauge action and the unimproved Wilson quark action (L « 3 fm). They 
reported that their data are described by these formulae. However, their sea quark masses are not quite small, and, 
since large scaling violation is susp ected with unimproved actions at coarse lattice spacings and lattice artifacts are 
suggested at strong couplings |23j, it should be demonstrated at weaker couplings in order that the discretization 
effects are actually under control. The UKQCD collaboration reported a result at mps/my = 0.44(2) obtained with 
the actions and the lattice spacing the same as those of JLQCD, with L w 1.6 fm |24|. They indicated the pion decay 
constant to bend slightly downward at this quark mass, but further work is required for quantitative comparison with 
the ChPT predictions. 

In this paper, we follow up on our previous two-flavor full QCD work 3] with an RG-improved gauge action and 
tadpole-improved O(o)-improved Wilson-clover quark action at mps/my — 0.80-0.55 and attempt to lower the quark 
mass to give mps/my down to 0.35. Since the computational costs grows rapidly toward the chiral limit, roughly 
proportional to (mps /mv)~ e we concentrate our effort on the coarsest lattice of a w 0.2 fm at = 1.8, while 
using improved actions. 

Generation of configurations below mps/my ~ 0.5 demands technical improvements. The BiCGStab algorithm 
sometimes fails to converge, which we overcome by an improvement called BiCGStab(DS-L) |2y,|27|. Another problem 
is the emergence of instabilities in the HMC molecular dynamics evolution [281 |29| . This seems to be caused by 
very small eigenvalues of the Dirac operator, leading to the change of the molecular dynamics orbit from elliptic 
to hyperbolic. The only resolution at present is to reduce the time step size. In this manner, we generated 4000 
trajectories at mps/my w 0.6, 0.5 and 0.4 and 1400 trajectories at the smallest quark mass of mps/my w 0.35 
on a 12 3 x 24 lattice with L w 2.4 fm. To examine the finite-size effect, we also generated 2000 trajectories at 
mps/mv ~ 0.6 and 0.5 on a 16 3 x 24 lattice with L w 3.2 fm. 

We calculate the light hadron spectrum and the quark mass on these configurations, and examine the validity of 
the quadratic chiral extrapolations by comparing the extrapolations made in the previous work with our new data at 
smaller quark masses. It turns out that the new data are increasingly lower than the extrapolation toward a smaller 
sea quark mass. We then examine how our data compare with the WChPT formulae, and whether WChPT fits using 
only the previous data at large quark masses predict correctly the new small quark mass data. This serves as a test 
to verify the viability of WChPT and of chiral extrapolations. 

Computing for the present work was made on the VPP5000/80 at the Information Processing Center of University 
of Tsukuba. We used 4 or 8 nodes, each node having the peak speed of 9.6 Gflops. The present simulation costed 
0.119 Tflops-years of computing time measured in terms of the peak speed. 

This paper is organized as follows. We describe configuration generations in Sec. [H] The method of measurement 
of hadron masses, decay constants, quark masses and the static quark potential is explained in Sec. IIIII The finite- 



3 



size effects on hadron masses are also discussed in the same section. Sec. IfVI discusses chiral extrapolations with 
conventional polynomials, and those based on ChPT are presented in Sec. [V] Our conclusion is given in Sec. I VII 
Preliminary results of these calculations were reported in Ref. j^lj ■ 



II. SIMULATION 



For the gauge part we employ the RG improved action defined by 

s <> = § l c ° E Ku l i x ) + c i E w™{x)\ . (i) 



The coefficients c = 3.648 of the lxl Wilson loop and c\ — —0.331 of the 1x2 Wilson loop are determined by an 
approximate renormalization group analysis |3l|. They satisfy the normalization condition cq + Sci — 1, and (3 = 6/g 2 . 
For the quark part we use the clover quark action |32j | defined by 

S q = X! 7 / ' P ' (2) 

D x , v = 5 xy - {(1 - 'Y fl )U x ,ft.8x+(i,y + (1 + 7fi)Ul, u dx,v+fi} ~ <W c swft y~] ^F^, (3) 

where k is the hopping parameter, the standard clover-shaped lattice discretization of the field strength and 
dpj, = (*/2)[7|t)7i/]- F° r the clover coefficient we adopt a meanfield improved value csw = u q 3 E3 where 

u = (W lxl ) 1/4 = (1 - 0.8412/T 1 ) 174 , (4) 

using the plaquette VF lxl calculated in one- loop perturbation theory [3l]]. This choice is based on our observation 
that the one-loop calculation reproduces the measured values well [34j. 

Our simulation is performed at a single value of j3 = 1.8 using two lattice sizes 12 3 x 24 and 16 3 x 24 to study finite 
size effects. The lattice spacing fixed from m p at the physical sea quark mass is 0.2 fm. We adopt four values of the 
sea quark mass corresponding to the hopping parameter n sea = 0.14585, 0.14660, 0.14705 and 0.14720. This choice 
covers mps/my = 0.60-0.35, extending the four values n sea = 0.1409, 0.1430, 0.1445, and 0.1464 corresponding to 
mps/mv = 0.80-0.55 studied in Ref. The simulation parameters are summarized in Tabled where we also list 
the number of nodes (PE's) employed and the CPU time per trajectory. Gauge configurations are generated using 
the Hybrid Monte Carlo (HMC) algorithm [35l l3r|. The trajectory length in each HMC step is fixed to unity. We 
use the leap-frog integration scheme for the molecular dynamics equation. 

The even/odd preconditioned BiCGStab [3]} is one of the most optimized algorithms for the Wilson quark matrix 
inversion to solve the equation D xy G y = B x . However, BiCGStab sometimes fails to converge at small sea quark 
masses. While the CG algorithm is guaranteed to converge, it is time-consuming. We find that the BiCGStab(L) 
algorithm |38j , which is an extension of BiCGStab to L-th order minimal residual polynomials, is more stable |27| . 
Figure ^ illustrates for a very light valence quark mass corresponding to mps/my = 0.27 that the BiCGStab(L), 
while not convergent for L = 1 and 2, succeeds to find the solution for L = 4. In practice, however, too large L also 
frequently introduces another instability from possible loss of conjugacy among the L vectors. The optimum value of L 
depends on simulation parameters. To avoid a tuning of L at each simulation point, we employ the BiCGStab(DS-L) 
algorithm 26] . This is a modified BiCGStab(L) in which a candidate of the optimum L is dynamically selected. We 
find that BiCGStab(DS-L) is much more robust than the original BiCGStab at small quark masses. We also find that, 
at large quark masses where the conventional BiCGStab converges, the computer time required for BiCGStab(DS-L) 
is comparable. See Fig. [21 Therefore, we adopt BiCGStab(DS-L) at all values of our sea quark masses. 

We employ the stopping condition | \DG — B\ \ < A in HMC. The value of A in the evaluation of the fermionic force 
is chosen so that the reversibility over unit length is satisfied to a relative precision of order 10 -8 or smaller for the 
Hamiltonian, 

|A.0| = \H reverse d — Hq\, (5) 

where H reverse d is the value of the Hamiltonian obtained by integrating to t = 1 and integrating back to t = 0. We 
also check the reversibility violation in the link variable, 




^l = J^ U;™p ed (n) - U°^ b (n), (6) 
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where the sum is taken over all sites n, colors a, b and the link directions fi. We illustrate our check in Fig. 
where results at n sea — 0.14585 and n sea = 0.14705 on 20 thermalized configurations separated by 100 trajectories 
are shown. When the sea quark mass is large (n sea — 0.14585, mps/mv — 0.6), the violation does not show any 
clear dependence on the stopping condition. For small sea quark mass (n sea = 0.14705, mps/my — 0.4), however, 
it depends on the stopping condition significantly. We must be careful with the choice of the stopping condition at 
small sea quark mass. We use a stricter stopping condition in the calculation of the Hamiltonian in the Metropolis 
accept/reject test. Table [I] shows our choice of A together with the average number, Ni nv , of the BiCGStab(DS-i) 
iterations in the quark matrix inversion for the force calculation. 

In the course of configuration generation by the HMC algorithm, we sometimes encountered extremely large values 
of dH = H t r i ai — Hn, the difference of the trial and starting Hamiltonians. Similar experiences have been reported by 
other groups [28ll29j| . Empirically this phenomenon occurs more frequently for smaller sea quark masses at a fixed step 
size, and can be suppressed by decreasing the step size. A typical example is shown in Fig. ^ In our runs we employ 
a step size dt small enough for this purpose. As a consequence our runs have a rather high acceptance 80-90%. 
It is possible that this phenomenon is connected to the appearance of very small eigenvalues of the Wilson-clover 
operator toward small quark masses. In the right panel of Fig.^J we show the norm ||.D~ 1 (Z^)~ 1 (/>|| (triangles) and 
the contribution of the smallest eigenvalue of 75D to the norm (filled squares). We observe that the jump of dH 
(open circles) is associated with a peak of the norm, and that the peak is saturated by the contribution of the smallest 
eigenvalue. We suspect that such small eigenvalues cause some modes of the HMC molecular dynamics evolution to 
change its character from elliptic to hyperbolic, leading to divergence of the Hamiltonian. We defer a further study 
of this problem to future publications. 

We accumulate 4000 HMC trajectories at n sea — 0.14585, 0.14660 and 0.14705 and 1400 trajectories at K sea = 
0.14720 on the 12 3 x 24 lattice. We also accumulate 2000 trajectories at n sea = 0.14585 and 0.14660 on the 16 3 x 24 
lattice. Measurements of light hadron masses and the static quark potential are carried out at every 5 trajectories. 



III. MEASUREMENT 
A. Hadron masses 

The meson operators are defined by 

M(x) = qM(x)Tqto(x), T = I, 7s,7„,757m , (?) 
where / and g are flavor indices and x is the coordinates on the lattice. The octet baryon operator is defined as 

O f9h (x) = e abc (V /)Q 0E) T CW ff)& (z)) q {h)c {x), (8) 

where a, b, c are color indices and C = 7472 is the charge conjugation matrix. Decuplet baryon correlators are 
calculated using an operator defined by 

Dl gh (x) = e abc ^^(xfCj^ix)) q (h)c (x). (9) 

For each configuration quark propagators are calculated with a point and a smeared source. For the smeared source, 
we fix the gauge configuration to the Coulomb gauge and use an exponential smearing function ip(r) = Aexp(-Br) 
for r > with ^(0) = 1. We chose A — 1.25 and B — 0.50 as in our previous study 3]. In order to reduce the 
statistical fluctuation of hadron correlators, we repeat the measurement for two choices of the location of the hadron 
source, t src — 1 and N t /2 + 1(= 13) and take the average over the two |ll| : 

\ [(H(t src + t)H(t src ^) tm=1 + (H(t src + t)H(t src )l) tarc=Nt/2+1 ) . (10) 

This procedure reduces the statistical error of hadron correlators typically by 30 to 40%, which suggests that the 
statistics are increased effectively by a factor of 1.7 to 2. For further reduction of the statistical fluctuation, we take 
the average over three polarization states for vector mesons, two spin states for octet baryons and four spin states for 
decuplet baryons. 

Figures |S] and illustrate the quality of effective mass plots. For mesons, an acceptable plateau of the effective 
mass is obtained from hadron correlators with the point sink and the doubly smeared source. Signals are much worse 
for baryons. 
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We carry out x 2 fits to hadron correlators, taking account of correlations among different time slices. A single 
hyperbolic cosine form is assumed for mesons, and a single exponential form for baryons. We set the lower cut of 
the fitting range as t m ; m — 6 for mesons and t m i n = 5 for baryons, which is determined by inspecting stability of the 
resulting mass. The upper cut (t max ) dependence of the fit is small and, therefore, we fix t max to N t /2 for all hadrons. 
Our choice of fit ranges and the detailed results of hadron masses are given in tables of Appendix^] Statistical errors 
of hadron masses are estimated with the jack-knife procedure. We adopt the bin size of 100 trajectories from an 
analysis of the bin size dependence of errors as discussed below in Sec. IIII El 



B. Quark masses 



We calculate the mean up and down quark mass through both vector and axial-vector Ward identities. The two 
types of quark masses, denoted by mvwi and itiawi respectively, differ at finite lattice spacings because of explicit 
violation of chiral symmetry by the Wilson term. 

A bare VWI quark mass is defined by 



1 fl 1 



K, 



The critical hopping parameter k c is determined by chiral extrapolations as discussed in Sec. HVI andlVl A bare AWI 
quark mass is calculated using the fourth component of the improved axial-vector current 

4 mp = A 4 + c A d 4 P, (12) 

where P is the pseudoscalar meson operator, Eq. Q with T = 75, and 84 is the symmetric lattice derivative. Then, 



rn 



AWI 
guar 



k is obtained through 



awi _ m PsC% 

l quark ~ ^(js ' V °J 



The amplitudes C A and Cf> are calculated as follows. We determine the pseudoscalar meson mass mps and Cp by 

(P'(t)P s (0)l) = C P [exp(-TOp S <) + exp(-mps(L t - t))} , (14) 

where the superscripts I and s distinguish local and smeared operators. Keeping rnps fixed, we extract C A from 

(4 mp ''(*)P s (0)t) = C% [exp(~m PS t) - cxp(-m PS (L t - t))} . (15) 

The renormalized quark masses in the MS scheme at 2 GeV are obtained as follows. The VWI up and down quark 
mass 



ywi _ 1 ' 1 ' 



l U(l n I 



(16) 



with k U 4 the hopping parameter at the physical point, is renormalized using one- loop renormalization constants and 
improvement coefficients at fl = 1/a: 

<T IM % - V«) = Z m ( 1 + bj^- ) 2^L_. (17) 



U J Uq 



Similarly the renormalized AWI quark mass is obtained by 



— Z A [l + b A - 

AWIMS/ t / \ \ u o J AWI / 1C A 

m ud ' (/x - 1/a) = 7 m vwA m ud , (18) 

Zpil + h 



where is the value of m^^ k extrapolated to n u d- The determination of n u d is discussed in Sec. II VI and Ivl 

Since non-perturbative values for the renormalization coefficient Z A and the improvement parameters c A , b A etc. arc 
not available for our combination of actions in two-flavor QCD, we adopt one-loop perturbative values calculated in 
Refs. [23, |4(| improved with the tadpole procedure using uq given in Eq.QJ. The MS quark masses at // = 1/a are 
evolved to = 2 GeV using the four- loop beta function |4ll l42j . 
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C. Decay constants 

The pseudoscalar meson decay constant is calculated by 



/ m vwi \ rs r i 

f P8 = 2nu Z A 1 + bA ^± £4 J ±£_ (19) 
^ uq J Gf, y mps 

where C l P is determined by 

(P*(t)P'(0)t) = C l P [cM-m P st) + exp(-m PS (i t - *))] , (20) 

keeping mps fixed to the value from (P l (t)P s (0) <! ). 
The vector meson decay constant fy is defined as 

(0\Vi\V) = eifvmv, (21) 

where is a polarization vector. The procedure to obtain the vector meson decay constant is parallel to that for fps- 
The vector meson correlator with a smeared source is fitted with 

(V l (t)V s (0)^) = C v [exp(-mvt) + exp(-m y (i t - t))} , (22) 

which determines my and Cy. Using my as an input we fit the correlator 

{V l (t)V l (Q)1) - C l y [exp(-myt) + exp(-m y (i 4 - i))] , (23) 

where the amplitude C l v is the only fit parameter. A renormalized vector meson decay constant is then obtained 
through 



f v = 2 K u Z v [l + b v ^^-)J^, (24) 



u o / V m v 



where we also use one-loop perturbative values for Zy and by |39j, |40j . We do not include the improvement term 
Cvd v T nflu because the corresponding correlator is not measured. 

D. Static quark potential 

We calculate the static quark potential V(r) from the temporal Wilson loops W(r, t) 

W(r,t) = C(r)exp(-V(r)t). (25) 

We apply the smearing procedure of Ref. |43j ■ The number of smearing steps is fixed to its optimum value N opt — 2 
at which the overlap to the ground state C(r) takes the largest value. Let us define an effective potential 

V eff (r, t) = log [W(r, t)/W(r, t + 1)]. (26) 

Examples of V e ff are plotted in Fig. from which we take the lower cut of t min = 2. As shown in Fig. [SI we do not 
observe any clear indication of the string breaking. Therefore, we carry out a correlated fit to V(r) = V e ff(r,t m i n ) 
with 

V(r) =V -- + ar. (27) 
r 

Here we do not include the lattice correction to the Coulomb term calculated perturbatively from one lattice-gluon 
exchange diagram [44 , since rotational symmetry is well restored for our RG-improved action. The Sommer scale ro 
is defined through 45] 

4 ^ = L65. ,28) 
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We determine rg from the parameterization of the potential V(r): 



/ 1.65 — a ,„„. 

110 = V— 5— (29) 

The lower cut of the fit range in Eq. 1)2 7|l is determined as r m i n — \J~2 from inspection of the r m i n dependence of 
r . With 

Train ^ v2? / dof takes an unacceptably large value, while <x becomes ill-deterniined with T m i n > V~3. On 
the other hand, the r max dependence of ro is mild. Therefore, we fix r max to N s /2. We estimate the systematic error 
of the fit as follows. The fit of Eq. 1)27(1 is repeated with other choices of the range: t m i n = 3 or r m in = V3- The 
variations in the resulting parameters and ro are taken as systematic errors. The parameters in Eq. (|27f) and ro are 
presented in Table 1TT1 



E. Autocorrelation 



The autocorrelation in our data is studied by the cumulative autocorrelation time 

rg im (At max ) = - + £ Po(At), (30) 

At=l 

where po(t) is the autocorrelation function 

Po{At) = Y^' r o (At) = {(G(t)-{0))(0(t + At)-{0})}. (31) 

A conventional choice for At max is the first point where po vanishes because po should be positive when the statistics 
are sufficiently high. We take At max — 50 from the plaquette shown in Fig. |5J In Table IIIII we give T§ im for (i) 
the plaquette which is measured at every trajectory, (ii) the pseudoscalar meson propagators at t = N t /A, and (iii) 
the temporal Wilson loop with (r, t) — (2,2). Fig. 1101 shows the autocorrelation time for the plaquette. Combining 
the previous (open circles) and the new (filled circles) data, we observe a trend of increase for smaller quark masses. 
A sharp rise expected toward the chiral limit, however, is not seen. Our statistics may not be sufficient to estimate 
autocorrelation times reliably near the chiral limit. 

The bin size dependence of the jack-knife errors of hadron masses and Wilson loops is exhibited in Figs. ITTl The 
jack-knife errors reach plateaus at bin size of 50-100 trajectories. The situation is similar on 16 3 x 24. Therefore, we 
take the bin size of 100 trajectories in the error analysis. 



F. Finite-size effects 



In Figs. El an d El we present meson and AWI quark masses as a function of the spatial volume. The results 
obtained on 12 3 x 24 and 16 3 x 24 lattices are mutually consistent within errors. For baryons, there may be some 
indication in our data at mps/my = 0.50 (n sea — 0.14660) that the light baryon masses mjv and tua decrease by 
1-3% (0.8-3. la) as shown in Fig. El The effect is only around 2a, and higher statistics are needed to confirm if the 
difference can be attributed to finite-size effects. Finite-size effects in ro are expected to be much smaller than those 
in hadron masses. Our results in Fig. 1 151 confirm this. In the following analysis, we use data obtained on the 12 3 x 24 
lattice. 



IV. CHIRAL EXTRAPOLATION WITH POLYNOMIALS 



Extrapolation of the lattice simulation data to physical values requires some parameterization of the data as 
functions of the quark mass. In this section, we employ polynomials in quark masses. We work with the two data 
sets, the one obtained in the previous work that covers mps/my — 0.80-0.55 (the large quark mass data set), and 
the other obtained in the present work that covers mps/my = 0.60-0.35 (the small quark mass data set), and with 
the combined data set of the two. For the large mass data set we borrow the fit from the previous work. 

We fit hadron masses in lattice units rather than those normalized by ro. With our choice of the improved actions, 
ro exhibits only a mild sea quark mass dependence as shown below in Sec. II V CI and hence introducing ro does not 
change convergence of chiral extrapolations. From practical side, ro suffers from a large systematic error on coarse 
lattices with a = 0.2 fm. Hence fits become less constraining if hadron masses are normalized by ro- 
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A. Pseudoscalar meson mass and AWI quark mass 

A quadratic form fitted well our previous lattice data of the pseudoscalar meson mass with a reasonable x 2 /dof ~ 
1 Q. As shown in Fig. 1161 however, our new data at small sea quark masses deviate significantly from the quadratic 
fit. Inclusion of the small quark mass data set in the quadratic fit rapidly increases \ 2 /dof to ~ 10. In addition, the 
determination of the critical hopping parameter k c becomes unstable as shown in Fig. 1171 A reasonable /dof and 
a stable fit are achieved only when we extend the polynomial to quartic, 

2 — rPS^VWI , r PS/VWI \2 , n PS(VWI \3 , pPS/VWI \4 /on\ 

m PS - B m quark + C {m quark ) +D {m quark ) +h (m quark ) , (32) 

where rn^J^rk ^ s given in Ea. (|ll[l and n c is taken as a fit parameter. The quartic polynomial provides the best fit 
among our tests varying the order of polynomials. 

Since m 2 PS may be affected by the logarithmic singularity of ChPT, we examine the convergence of extrapolations, 
i.e., whether it depends on the order of polynomials, using Tn q ^^ k that has no logarithmic singularities. Along with 
the case of m 2 PS , the new data at small quark masses deviate from the quadratic fit obtained from the large quark 
mass data, as depicted in Fig. ^] We fit m q ^ k by 

^AWI _ nAWIVWI , r AWI (VWI \2 , r>AWI (VW1 \3 , pAW I / VW I \4 /oo\ 
171 quark — a m quark ' ° \ m quark) "+" u \ m quark) "+" & \ m quark) ■ 

The fit range and order dependence are given in Fig. 1191 (wij^4) 4 terms are needed again to obtain a reasonable 
X 2 /dof. 

We find that n c determined from mp s agrees with that from m q ^ k within errors. Hence we simultaneously fit 
m 2 PS and »nj^f fc to determine k c . The resulting independent and simultaneous fits to m PS and m q ^ k are presented 
in Tables llVl and Ivl respectively. The difference in mass from the fits including (m^^ k ) 5 is taken as systematic 
errors. These errors represent only uncertainties within polynomial extrapolations. As shown in Sec. IV Bl WChPT 
fits sometimes lead to values beyond these systematic errors. 



B. Vector meson mass 



We fit vector meson mass with a cubic polynomial in mp s , 

m v =A V + B v m PS + D v m% s + F v m PS (34) 

with the results shown in Fig. [^H] and Table IVII As in the case of m 2 PS and m^^^ k , systematic deviations from the 
previous fit are observed, although the difference (7% or 3.6<r in the chiral limit) is smaller. Inclusion of terms m PS 
and rn PS gives a good fit with a satisfactory Q. We estimate the systematic error from higher order terms by the 
difference from the fit with m PS term. 

The effects of vector meson decays are not considered in the fit. If a vector meson decays into two pseudoscalar 
mesons, a vector meson with the momentum p = 2tt/L will take a different energy depending on whether it is 
polarized parallel or perpendicular to the momentum direction, because of mixing of one vector meson state and two 
pseudoscalar meson state 0, 01 • We find no indication of vector meson decays as shown in Fig. |2] Our sea quark 
masses and the lattice size do not seem to be enough to allow the decay. 



C. Decay constants, baryon masses and Sommer scale 

Chiral extrapolations are carried out for pseudoscalar and vector meson decay constants and octet and decuplet 
baryon masses using cubic polynomials in rn ps , 

f PSV = A fpsJv + B fpsJv m PS + D fps - fv m PS + F fpsJv m PS , (35) 
m oct , d ec = A oct ' dec + B oct ' dec m 2 PS + D oct ' dec m PS + F oct ' dec m% s . (36) 

The results are presented in Figs. 1221 and!23l and Tables IVTTl and fVIIII While the decay constants show clear deviations 
from the previous fit, baryon masses are almost on the fit. We gather that the latter is an accidental effect that is 
caused by a compensation of the downward shift of baryon masses expected toward a small quark mass with an 
upward finite-size shift caused by somewhat too small a lattice (L = 2.4 fm) for baryons (see Sec. IIII F]l . 
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The Sommer scale tq is often extrapolated linearly in ra 2 PS . Since we find a curvature in our data, however, we 
adopt the same form as that for the vector meson masses, 

— = A ra + B r °m 2 PS + D r «m PS + F r °m PS . (37) 

The results are seen in Fig. |2] and Table |I 

D. Results at the physical point 

The physical point is defined by empirical pion and p meson masses, M r = 0.1350 GeV and M p = 0.7711 GeV. 
With our polynomial fit, the physical point m n for mps is determined by solving the equation, 

M w (3g) 



A v + B v ml + D v m% + F v m% M p 



The p meson mass at the physical point m p is obtained by Eq. 1341) with nips = m^, which determines the lattice 
spacing a p = 0.2007(38) fm. The lattice spacing can also be determined from ro taking its phenomenological value 
R = 0.49 fm. Using Eq. instead of Eq. we have 

= M«R . (39) 



Substitution of to Eq. I|37|l leads to ro at the physical point, yielding an alternative lattice spacing a ro , a ro = 
0.2119(61) fm, which is consistent with a p within 2a. 

We calculate ra^ 1 using n u d defined by mps{K u d) — m^, and vn^ 1 by Eq. I|33l) . and then convert to renormalized 
quark masses in the MS scheme at 2 GeV (see Sec. IIIIBfl . Table |X| presents a summary of the parameters at the 
physical point, obtained with polynomial extrapolations, together with comparisons with the quadratic fit in the 
previous work. The difference between old and new results is generally 4-8% except for the VWI quark mass for 
which a difference more than 20% is observed (see Fig. 12511 . The latter is caused by a shift of k C) with which even a 
small shift leads to an amplified change in the mean up and down quark mass. 

V. CHIRAL EXTRAPOLATION BASED ON CHPT 

We first examine the one-loop formulae from continuum ChPT, which have already been tested in 0,|23|. We 
then attempt a fit based on WChPT including effects of 0(a 2 ) chiral symmetry violation due to the Wilson term. 

A. ChPT extrapolation 

The one-loop formulae |l2lll4j| derived from ChPT in the continuum limit are 

= 1+ lQ g Tl ( 40 ) 



2B m quark 2 (4^/) 2 6 A| 

„ » / 1 2Born quar k 2Bom quar k\ .. . 

fps - f [ 1 ~~(Mf) S ~ S A! J ! (41) 

where Bq, f, A3 and A4 are parameters to be obtained by fits. The coefficient 1/2 in front of the logarithm is 
a distinctive prediction of ChPT. Since several parameters are common in the two formulae, we fit m 2 PS and fps 
simultaneously. Correlations between m 2 PS and fps are neglected in the fits for simplicity. Thus, the x 2 /dof serves 
only as a guide to judge the relative quality of the fits. We estimate the errors by the jackknife method. We try both 
m quark an d m JLT/fc (Cases 1 and 2 in what follows) for m quar k that appears in these formulae. For m quar k = m q / mrki 

AWI 
quark 

in Table KH we find: 



we use k c determined in Eq. I|33l) since rn^^^ k has no logarithmic singularities in ChPT. From the fits summarized 



Case 1 {m qU ark = m q ^rk) '■ When we fit the data over the whole range mps/mv = 0.80-0.35, we are led to a 
large y^/dof ~ 70. By restricting the fitting interval to mpsjmy — 0.60-0.35 we obtain a reasonable fit 
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with x 2 /dof — 1.9, which is plotted in Fig. [25] As one observes in the second panel of this figure, which 
shows m%s/^ m ^uark appearing in the left hand side of Ea. l|4U|l . the chiral logarithm may be visible only at 
mps/mv < 0.40. 

Case 2 (m quar k — m q^rk) '■ ^ n contrast to Case 1, m ps /2mY^ k increases toward the chiral limit in the whole mass 
range, which is seen in Fig.|57| Nevertheless, the situation is similar. A fit over the whole range mps/mv = 0.80- 
0.35 leads to x 2 /dof ~ 100. To obtain an acceptable fit, we have to remove the data at large quark masses. 
The best fit obtained for the range mps/mv = 0.60-0.35 is shown in Fig. 1271 

In neither case do we draw the clear evidence for the chiral logarithm for pseudoscalar mesons. 
For the vector meson, we adopt the formula based on ChPT in the static limit . 



my 



A v + B v m 2 PS + C v m PS . (42) 



This cubic form describes our data well as shown in Fig. [281 (see Table Kill for numbers). 
For octet and decuplet baryons we employ a similar cubic formula .49] 



a oct, dec , r>oct,dec 2 , sioct,dec'i ( A*)\ 

m oc t,dec — A + u m PS + C m PS , [46) 



which also reproduces our data well (Fig. and Tabic 1X111(1 . 

In order to present predictions at the physical point, we carry out extrapolations using the data at mps/rnv = 0.60- 
0.35. From Eq. I|42|) the physical point for m P g is given by 

M - (44) 



AY + B v ml + C v ml M p 



The lattice spacing is determined to be a^ hPT — 0.192(10) fin. For the vector meson, a fit for the whole range 
m PS /m v = 0.80-0.35 is acceptable, as seen in Table ED We will use this fit in Sec. IVBl with a c p hPT = 0.2009(21) fin 
for this case. 

The masses of non-strange baryons N and A are determined by substituting m, to m P s in Eq. 143(1 . The bare quark 
mass at the physical point m u d and the pion decay constant f v are obtained from Eqs. (|4t)|l and 1(41(1 . Renormalized 
quark masses are calculated with m u d as in the case of polynomial extrapolations. These results are compiled in 
Table IXIVl 

We observe 5-10% differences between the ChPT fits over mpsjrny — 0.60-0.35 and the quadratic polynomial 
fits over mps/my = 0.80-0.55 obtained in the previous work. The numbers are tabulated in Table IXl These 
differences are similar in magnitude to those we found with higher order polynomial extrapolations using the whole 
range mps/my — 0.80-0.35. An exception is the VWI quark mass on which we shall make a further comment below. 



B. WChPT extrapolation 



WChPT without resummation 



ChPT adapted to Wilson-type quark actions on the lattice (WChPT) has been addressed in Refs. 0, 0, . 
An important point i s that 0(a 2 ) chiral breaking terms in the chiral Lagrangian are essential to generate the 
parity-flavor breaking phase transition , which is necessary to explain the existence of massless pions for Wilson- 
type quark actions |5ft l5lt IH^. Therefore, we must include the 0(a 2 ) terms in the leading order. In this counting 
scheme, the one-loop formulae read |l8|. 



mp s 



in 



AWI 

quark 



IPS 



Am VWI 1 4-,,, PS m VWI lotr 
/im quark 1 + u l m quark I0 6 



/ / Am VWI 

m VWI -I ,,.AWI VWI i "V^ 
m quark I 1 + w l m quark 10 S I ^2 

/ / Ain VWI 

f 1 -~iJ PS rn VWI Iop- 1 uark 
/I 1 w l m quark 10 g I ^2 




(45) 
(46) 
(47) 
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Here k c in m^]^/ fc , A, f, ljq, uf s , uj^ wi , uj[ f s , A3, A 3 ^awi, A and A 4 are free parameters, and the overall factor of 
m quark 1S absorbed in ujq and Ao. We note that A consists of O(a ) and 0(a 1 ) parts, ujq ~ 0(a 2 ), lo^ wi ~ O(a) and 

PS _ l (A + <o) 
Wl " 2l4^)^' (48) 

1 " (Wo) 2 ' m 

where fo is the pion decay constant in the continuum and chiral limit, which can be different from / by 0(a). The 
constants w± and are O(a ). 

There are two features in these formulae worth emphasizing. First, the coefficients of mq Ua rk^ogm qua rk terms 
receive contributions of 0(a). This is in contrast to continuum ChPT, in which these coefficients take universal 
values. Second, there are terms of the form a 2 \ogm quar k which are more singular than the m quar k^ogm quar k terms 
toward the chiral limit at a finite lattice spacing. Thus WChPT formulae predict the chiral behavior at finite lattice 
spacings that is different from what is expected from ChPT in the continuum limit. 

We fit mps and rn^^^ k simultaneously, neglecting correlations between them. The errors are estimated by the 
jackknife method. We then fit fps with A and k c fixed from Eqs. I|45(l and (|46|l . We give the results in Fig. 1301 and 
Tables IXVl and IXVfl Fig. QUI demonstrates that the one-loop WChPT formulae explain our data over the whole range 
mps/mv = 0.80-0.35. 



2. Resummed WChPT 



While fits with Eqs. (|45|l and i|46|) work well for the whole range of quark mass we measured, extrapolation to the 
physical point is still problematic because the ujq log rn^^rk terms become larger than the leading terms in the chiral 
limit. A way out has been proposed in Ref. fl8| in which leading singularities around the chiral limit is resummed. 
The resulting formulae read, 



l PS 



AWI 

quark 



Am VWI 

quark 



rn 



vwi 




ps^vwi 



rn 



qu.ar 



AWI VWI 



m vwl lnir 
m quark 10 S 



A 2 



A 2 

\ IY 3,AWI 



(50) 
(51) 



where the fitting parameters are k c in m^7 fc , 



A, ujq, u( b ,uj 



,AWI 



A3, A-3 t Awi and Aq. The minus sign in the resummed 



part is introduced to keep — log^m^^/A 2 ,) positive. We note that fps is not affected by the resummation except 
for a shift of k c . 

As with the case of WChPT without resummation, these resummed WChPT formulae describe our data for the 
whole range of mps/mv = 0.80-0.35. The results are seen in Fig. 1311 and Tables IXVIII and IXVIIII 

The magnitude of the leading and the one- loop contributions is plotted in Fig. 1321 as a function of m quar k ■ In 
contrast to WChPT without resummation, which is shown in the second panel of the figure, the one-loop contribution 
of resummed WChPT fit remains small in the the whole range of quark mass we explored, including the chiral 
limit. This confirms the convergence of the resummed WChPT formulae. Furthermore, the resulting parameters are 
comparable with phenomenological estimates; we obtain A3 = 0.15(15) [GeV]and A4 = 2.44(13) [GeV] as compared 
to A3 = 0.2-2.0 [GeV] and A4 = 1.26(14) [GeV], respectively, from Ref. [l2l |53| . A more accurate examination 
requires extrapolation to the continuum limit, which is left for studies in the future. 

In the present fit, the m qua rk ^ogm quar k terms are sizably suppressed due to O(a) corrections for the pseudoscalar 



K 



PS 



,AWI 



) represents the strength of the chiral logarithm. The 



meson mass. In the combination mp s / ^m^/j., 
resummed WChPT fit gives (wf s - lj^ wi ) = 0.24(13), while in continuum ChPT we expect s = A/32ir 2 f$ = 2.7 
and u!^ n 1 — 0, with the phenomenological value of fo — 0.086 GeV, ignoring 0(a) dependence in A. Namely, 
the coefficient of the logarithm is suppressed to about 10% of the ChPT value by 0(a) contributions in m PS and 
It is important to repeat a similar analysis at a smaller lattice spacing to verify that the magnitude of the 



AWI 
quark ' 



iquark log m quar k coefficient converges toward the value predicted by ChPT. 
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Results at the 



point 



Since WChPT formulae are not available for the vector meson, we adopt Eq. I|42[) to fix the physical point for . A 
fit for the whole data in the range mps /m v = 0.80-0.35 yields a c p hPT = 0.2009(21) fm. Substituting m T to Eq. (5QJ 
and using a^ hPT , we obtain the VWI quark mass at the physical point 



Eqs. O and O with 



,vwi 



then 

yield and /„■ respectively fTable |XIX)| . 

Let us compare the resummed WChPT results with those of the quadratic polynomial obtained with the original 
data over the range mps/my — 0.80-0.55 (Table EJ and the fits using ChPT formula in the continuum limit for 
mpsl m v = 0.60-0.35 fTable |XIV|) . The lattice spacing, the AWI quark mass, and the pion decay constant take 
similar values among higher order polynomials, ChPT and resummed WChPT formulae. An exception is the VWI 
quark mass which significantly depends on the functional forms for the chiral extrapolation, (see Fig. I33J1 . Our final 
values for the light quark mass at a = 0.2 fm are: 



VWI, MS I r, ^ Tr\ 

m„,,, \ji — 2 CeV) = 



in 



"ud 

_AWI,MS 
ud 



(p = 2 GeV) 



1.314(99) [MeV] 
1.796(51) [MeV] 

2.902(36) [MeV] 
2.927(53) [MeV] 



(resummed WChPT) 
(polynomial) 

(resummed WChPT) 
(polynomial) 



(52) 
(53) 



The sensitivity of the VWI quark mass on the functional form of chiral extrapolation is due to closeness of K u d to the 
critical value k c . A small variation of k c is easily amplified in the up and down quark mass which is determined by 
the difference 1/ ' k u( i — 1/k c . 



4- Chiral extrapolation from large quark masses 



Finally, we test if WChPT explains the deviations of our new data at small quark masses from the quadratic 
extrapolation of the data at mps/my — 0.80-0.55. A motivation of this test is the rapid increase of the computational 
time to simulate QCD toward small sea quark masses on fine lattices. If WChPT correctly predicts the small quark 
mass behavior from heavy sea quark mass simulations for mps/my > 0.5, it will be a great help for our studies. 

We apply the resummed WChPT formulae to the large quark mass data set at j3 = 1.80. Since the number of data 
points at m PS /m v > 0.5 is small for a stable fitting, we introduce a restriction: A 3 = A 3 awi- Fig.|^(see Table IXIXl 
for numerical values) compares the fit from the large quark mass data set and that using the data for the entire mass 
range. The resummed WChPT fit using the large quark mass data set alone describes the small sea quark mass data 
very well. This contrasts to the polynomial extrapolation. Our observation suggests that WChPT may provide a 
valuable tool to carry out an accurate chiral extrapolation using simulations with not too small quark masses. 

Encouraged by this, we apply the resummed WChPT to the two additional data sets at mps/my = 0.80-0.55 
obtained at smaller lattice spacings at p = 1.95 and 2.1 (a = 0.16 and 0.11 fm) in the previous work. A simultaneous 

linear continuum extrapolation using nz^' / ' MS and ni^^ 1 '^ , combined with the results for P — 1.8, leads to 

mjjjf (/* = 2 GeV) = 3.06(18) [MeV] (resummed WChPT fit) (54) 
where the error is statistical only. When we use our whole data of mps/my — 0.80-0.35 at P = 1.80, we obtain 

m^f {fi = 2 GeV) = 3.11(17) [MeV] (resummed WChPT fit with our whole data) (55) 

This is compared to our previous result using the quadratic extrapolation: 

m^f (jx = 2 GeV) = 3.45(10) [MeV] (quadratic fit) (56) 

The resummed WChPT results in a 10% decrease in the mean up and down quark mass. This is demonstrated in 
Fig. ESI 



VI. CONCLUSIONS 



In this paper, we have pushed our previous study of two-flavor QCD down to a sea quark mass as small as 
rnps/my = 0.35, using the RG-improved gauge action and the clover-improved Wilson quark action. We have 
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found that our new data at mpsjmy — 0.60-0.35 show clear deviations from the prediction of the previous chiral 
extrapolations based on quadratic polynomials, which implies that higher order terms were needed to describe the 
behavior at a small sea quark mass. On the other hand, our current data do not show the clear quark mass dependence 
expected from ChPT in the continuum: the chiral logarithm may appear only below mps/my ~ 0.4. This result 
contrasts with that of the qq+q collaboration 22] based on unimproved plaquette glue and Wilson quark actions, but 
is not dissimilar to that of UKQCD Q • 

We have provisionally ascribed the major reason for the failure of continuum ChPT to explicit chiral symmetry 
breaking of the Wilson term, which is significant on our lattice of a = 0.2 fm. We have then made a test of WChPT in 
which the effect of the Wilson term is accommodated, and found the resummed one-loop WChPT formulae that take 
account of the effects up to 0(a 2 ) describe well our entire data. Convergence tests indicate that resummed WChPT 
gives well-controlled chiral extrapolations. The use of WChPT generally leads to modifications of various physical 
observables at the physical point by about 10%, compared with those obtained in the quadratic extrapolation at this 
lattice spacing. A much larger modification, however, is seen with the light quark mass defined through vector Ward 
identity: the WChPT extrapolation decreases it by 30%. 

We note in particular that the resummed WChPT extrapolation from our previous data at mps/mv = 0.80- 
0.55 predicts correctly the new data at mps/mv = 0.60-0.35. Encouraged by this fact, we attempted a continuum 
extrapolation of the light quark mass using the resummed WChPT fits to the previous data at raps/my — 0.80-0.55 
but on finer lattices with a = 0.16 and 0.11 fm. We find in the continuum limit, m^ s (/i = 2 GeV) = 3.11(17) [MeV], 
which is smaller than the previously reported result by approximately 10%. Our work suggests that WChPT provides 
us with a valuable theoretical framework for chiral extrapolations. 
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APPENDIX A: HADRON MASSES 

Measured hadron masses are summarized in Tables IXXlHXX VlTl Our choice of the fitting range and resulting value 
of x 2 /dof are also given in these tables. 
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TABLE I: Run parameters of the present simulation. The step size dt is given by the inverse of the number of the molecular 
dynamics steps (#MD), and hence not listed. We denote the tolerance parameter in the stopping condition for the quark 
matrix inversion in calculations of the force by A/ orce and the average number of iterations by Ni nv . Number of node (PE's) 
of VPP5000/80 used for the present calculation, and the CPU time required per trajectory in units of hour are also given. The 
number of trajectory is denoted by N tra j- 



12 3 x 24 


on 4PE 






on 4PE 




on 4PE 




on 8PE 




0.14585 






0.14660 




0.14705 




0.14720 


#MD 


200 




333 




400 


800 


1000 


1250 1600 


Accept. 


0.76 




0.72 




0.84 


0.82 


0.90 


0.87 0.91 


Ntraj 


4000 




1750 




2250 


680 


3320 


100 1300 


\ J'orCG 


1Q -io 






io- 11 




io- 11 




io- 12 


N inv 


87 






147 




232 




318 


Hour/traj 


0.23 




0.56 




0.69 


2.0 


2.6 


2.2 3.2 


mps/mv 


0.609(2) 






0.509(5) 




0.413(8) 




0.349(19) 


16 3 x 24 


on 4PE 






on 8PE 












0.14585 






0.14660 




0.14705 




0.14720 


#MD 


200 


250 


333 




500 








Accept. 


0.61 


0.71 


0.79 




0.80 








Ntraj 


800 


1200 


325 




1675 






— 


^force 


1Q -io 






io- 11 










Ninv 


92 






158 










Hour/traj 


0.50 


0.61 


0.69 




1.03 








mps/mv 


0.604(3) 






0.509(4) 










TABLE II: 


String tension a and Sommer scale ro at simulated 


sea quark masses. The first c 


3rror is 


statistical. The second and 


third ones 


are the systematic errors due to the choice of tmin and r min . 








12 3 x 24 




train 




f*min 




a 






ro 


0.14585 


2 




V2 




0.322(6)(- 


-42) (+91) 




2.004(8)(+58)(+77) 


0.14660 


2 




V2 




0.289(5)(- 


•8) (+64) 




2.107(8)(+37)(+54) 


0.14705 


2 




V2 




0.278(5)(- 


■38) (+34) 




2.167(9)(+80)(+25) 


0.14720 


2 




V2 




0.255(8)(- 


-10)(+42) 




2.237(17) (+10) (+34) 


16 3 x 24 


f^sea 


tmin 




Tmin 




a 






ro 


0.14585 


2 








0.313(11)( 


-10)(+90) 




2.011(10)(+17)(+72) 


0.14660 


2 




^/2 




0.270(6) (+7) (+66) 




2.131(ll)(+20)(+39) 
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TABLE III: Autocorrelation time for plaquette (t^™), pseudoscalar meson propagator at N t /4 (r^g n ) and Wilson loop with 
(r,t) = (2,2) (r^ m ). All values are in units of HMC trajectory. 



12 3 x 24 


f^sea 


0.14585 


0.14660 


0.14705 


0.14720 


cum 
'plaq 


7.6(1.8) 


11.7(2.3) 


9.5(2.1) 


8.9(3.2) 


cum 
T PS 


7.9(1.6) 


7.2(1.5) 


5.3(1.2) 


3.0(1.0) 


cum 
T W 


8.1(1.9) 


12.6(2.9) 


11.3(2.2) 


13.0(4.4) 


16 3 x 24 




0.14585 


0.14660 


0.14705 


0.14720 


cum 
'plaq 


14.1(3.9) 


8.8(2.1) 






cum 
TPS 


10.3(2.8) 


4.9(1.6) 






cum 
T W 


14.1(3.8) 


10.1(4.3) 







TABLE IV: 


Parameters of indepene 


lent polynomial 


chiral fits to 


AWI quark masses 


and pseudosca 


lar meson masses 


as a 


function of the VWI quark mass. 














nips/mv 




gAWI 




jjAWI 




X 2 /dof 


Q 


0.80-0.35 


0.147502(14) 


1.961(60) 


— 10.5(1.9) 


71(20) 


— 201(67) 


4.38/3 


0.22 


mps/m v 




b ps 


c ps 


d ps 


e ps 


X 2 /dof 


Q 


0.80-0.35 


0.147514(15) 


12.05(33) 


-55.7(9.0) 


359(89) 


-966(281) 


4.17/3 


0.24 



TABLE V: Parameters of simultaneous polynomial chiral fits to AWI quark masses and pseudoscalar meson masses as a 
function of the VWI quark mass. The first error is statistical and the second is a systematic one due to the higher order term 
for the chiral extrapolation. 



mps/mv k c 


B AWI qAWI 


jjAWI 


gAWI 


X 2 /dof 


0.80-0.35 0.147508(14) (+7) 


1.938(54) (-60) -9.8(1. 7) (+3.3) 


65(18)(-67) 


-181(60)(+541) 


8.89/7 




b ps c ps 


d ps 


e ps 


Q 



12.18(31)(-20) -58.9(8.6)(+7.9) 389(85)(-126) -1053(269)(+880) 0.26 



TABLE VI: Parameters of polynomial chiral fits to vector meson mass. The first error is statistical and the second is a 
systematic one due to the higher order term for the chiral extrapolation. 



nips/mv 


A v 


B v 


D v 


F v 


X 2 /dof 


Q 



0.80-0.35 0.770(15)(-6) 0.790(72)(+47) -0.304(97)(-116) 0.063(39)(+113) 1.10/4 0.89 



TABLE VII: Parameters of polynomial chiral fits to pseudoscalar and vector meson decay constants. 



mps/mv 


A fps 


B Jps 


D fps 


pfps 


X 2 /dof 


Q 


0.80-0.35 


0.1239(26)(-83) 


0.165(17)(+86) 


-0.076(27)(-264) 


0.018(12)(+296) 


17.2/4 


0.0018 


mps/mv 


A fv 


B fv 


D fv 


pfv 


X 2 /dof 


Q 



0.80-0.35 0.228(12)(-15) 0.265(59)(+127) -0.156(85)(-336) 0.039(35) (+340) 2.31/4 0.68 



TABLE VIII: Parameters of polynomial chiral fits to octet and decuplet baryon masses. 



mps/mv 


j^OCt 


■QOCt 




j-yOCt 


poet 


X 2 /dof 


Q 


0.80-0.35 


1.051(23)(+21) 


1.41(12)(- 


-19) 


-0.51(18)(+56) 


0.097(78)(-591) 


4.96/4 


0.29 


mps/mv 


A icc 


B acc 




D dec 


F icc 


X 2 /dof 


Q 


0.80-0.35 


1.381(37)(+8) 


1.03(19)(- 


-7) 


-0.17(29)(+18) 


-0.022(126)(-186) 


0.75/4 


0.95 
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TABLE IX: Parameters of polynomial chiral extrapolation of vq. 



mps/mv A rg 


B r(j D r() 


Fr X 2 /dof 


Q 


0.80-0.35 0.432(12)(- 


7) 0.207(89)(+99) -0.11(14)(-32) 


0.032(64)(+364) 0.37/4 


0.99 


TABLE X: Results of physical quantities obtained by polynomial chiral fits using data at mps/mv = 0.80-0.35. The results 
of the previous quadratic fits at mps/mv = 0.80-0.55 0] are also shown. The first error is statistical and the second is a 
systematic one due to the higher order term for the chiral extrapolation. Only statistical errors are given for the previous 
results. 




Quartic fit (this study) a 


Quadratic fit [3] 


Difference 


fit range in mps/mv 


0.80-0.35 


0.80-0.55 




a p [fm] 


0.2007(38) (-14) 


0.2150(22) 


-7%(6.5<j) 




0.147440(13)(+7) 


0.147540(16) 


0.1%(6.3cr) 


mr™(^2GeV)[MeV] 


1.796(51) (+18) 


2.277(27) 


-21%(18<r) 


m$Y lMa {.V> = 2 GeV)[MeV] 


2.927(53)(-55) 


3.094(35) 


-6%(4.8<r) 


MGeV] 


0.1248(31) (-59) 


0.1288(33) 


-3%(1.2<r) 


/p[GeV] 


0.2294(74)(-lll) 


0.2389(47) 


-4%(2.0cr) 


mjv[GeV] 


1.060(27) (+24) 


1.016(16) 


f4%(2.8a) 


m A [GeV] 


1.377(39) (+16) 


1.270(23) 


f8%(4.7a) 


"For vector meson masses, decay constants and baryon masses, we employ cubic fit functions in inpg as Eqs. 1341—1^(11. 




TABLE XI: Chiral extrapolation of pseudoscalar meson masses and decay constants based on the continuum ChPT formulae 
at one- loop with m quar -k = m^™' k and m quar k = rn^^rk- K c has been determined with mf^' k . The value of k c is shown in 
Table HVI 


mps/mv B 


pAWI j^AWI 


A AW I 2 I i f 

A 4 x /dof 


Q 


0.80-0.35 3.838(15) 
0.60-0.35 3.398(52) 


0.12162(47) 1.553(10) 
0.1130(20) 0.902(71) 


2.633(15) 849/12 
2.591(98) 11.2/6 


1Q -174 

0.083 


mps/mv B 


jVWI j^VWI 


A 4 ™ 7 xVdof 


Q 


0.80-0.35 6.886(22) 
0.60-0.35 6.582(87) 


0.13225(35) 2.4018(85) 
0.1145(18) 1.645(83) 


2.463(11) 1417/12 
2.262(72) 17.5/6 


1Q -296 

0.0076 



TABLE XII: Parameters of chiral fits to vector meson mass based on continuum ChPT. 



mps/mv 


A v 


B v 


C v 


X 2 /dof 


Q 


0.80-0.35 


0.7692(86) 


0.897(32) 


-0.346(23) 


1.39/5 


0.93 


0.60-0.35 


0.731(45) 


1.31(49) 


-0.85(60) 


0.33/2 


0.85 



TABLE XIII: Parameters of chiral fits to octet and decuplet baryon masses based on continuum ChPT. 



mps/mv 




B oct 




X 2 /dof 


Q 


0.80-0.35 


1.043(14) 


1.641(68) 


-0.632(51) 


5.13/5 


0.40 


0.60-0.35 


1.011(52) 


2.08(59) 


-1.23(74) 


2.23/2 


0.33 


mps/mv 




Qdec 




X 2 /dof 


Q 


0.80-0.35 


1.351(20) 


1.353(88) 


-0.481(66) 


1.24/5 


0.94 


0.60-0.35 


1.428(86) 


0.52(93) 


0.53(1.15) 


0.23/2 


0.89 



18 



TABLE XIV: Results of physical quantities obtained by continuum one-loop ChPT chiral fits using data at mps/mv = 0.60- 
0.35. For the m quar k = m^^'k case, k c has been fixed to the value determined from the quartic fit to mf^' k shown in Tablc lTvl 
The errors are statistical. 





continuum ChPT (m quark = mf^'k) 




continuum ChPT (m quark = m%™* h ) 


dp [fin] 


0.192(10) 




0.147445(14) 




0.1474431(65) 


VWIMSr 
m ud \M 


= 2 GeV)[MeV] 


1.609(89) 




1.625(81) 


AW I, MS/ 


= 2 GeV)[MeV] 


2.66(13) 




2.68(13) 


WGeV] 


0.1219(64) 




0.1231(65) 


m N [GeV] 


1.074(69) 


m A [GeV] 


1.47(11) 



TABLE XV: Parameters of chiral fits to pseudoscalar meson and AWI quark masses based on WChPT. 



mps/mv k c 


A 


LOO 


PS 


AWI 
U3 X 


Ao 


A 3 


Af WI 


X 2 /dof Q 



0.80-0.35 0.147445(27) 6.312(44) -0.40(13) -2.0(1.4) -2.0(1.4) 0.91(35) 1.95(15) 1.77(23) 11.9/8 0.16 



TABLE XVI: Parameters of chiral fits to pseudoscalar meson decay constants based on WChPT. n c and A has been fixed to 
the values in Table IXVl 

mps/my / u{ PS A4 X 2 /dof Q 

0.80-0.35 0.1233(17) 3.73(30) 2.44(13) 18.1/5 0.0028 



TABLE XVII: Parameters of chiral fits to pseudoscalar meson and AWI quark masses based on the resummed WChPT. 

mps/my k c A ua> lo? s u? WI A A 3 Aj wl x 2 /dof Q 

0.80-0.35 0.147459(20) 6.354(59) 0.542(46) 0.65(51) 0.42(49) 0.397(56) 0.15(15) 0.07(16) 11.0/8 0.20 



TABLE XVIII: Parameters of chiral fits to pseudoscalar meson decay constants based on the resummed WChPT. k c and A 
has been fixed to the values in Table IXVIII 



mps/my / lo{ ps A4 x /dof Q 

0.80-0.35 0.1227(17) 3.78(30) 2.44(13) 18.2/5 0.0028 



TABLE XIX: Results of physical quantities obtained by the resummed WChPT fits using data at mps/mv = 0.80-0.35. The 
results are compared with the results of the resummed WChPT fits using our previous data at mps/mv = 0.80-0.55. 





RWChPT 


RWChPT 11 


Difference 


fit range in mps/mv 


0.80-0.35 


0.80-0.55 




a P [fm] 


0.2009(21) 


0.2022(38) 


-l%(0.3a) 




0.147409(16) 


0.14736(22) 


+0.03%(0.2cr) 


™:r- MS (M = 2GeV)[MeV] 


1.314(99) 


1.10(64) 


+ 19%(0.3cr) 


™i™(M = 2GeV)[MeV] 


2.902(36) 


2.945(60) 


-1%(0.7<t) 


MGeV] 


0.1238(21) 


0.1368(43) 


-10%(3.0cr) 



"For mps/my = 0.80—0.55 data, we employ a restriction A3 = A^^awi- 
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TABLE XX: Meson masses and bare AWI quark masses on 12 3 x 24 lattice. 





mps 




X 2 /dof 


my 




X 2 /dof 


m AWI 

' ' L quark 


0.14585 


0.6336(14) 


[6,12] 


0.76(84) 


1.0405(38) 


[6,12] 


0.40(51) 


0.06340(34) 


0.14660 


0.4789(23) 


[6,12] 


1.60(1.19) 


0.9410(81) 


[6,12] 


2.36(1.02) 


0.03632(39) 


0.14705 


0.3520(29) 


[6,12] 


0.60(77) 


0.8526(148) 


[6,12] 


0.67(81) 


0.01952(30) 


0.14720 


0.2893(61) 


[6,12] 


0.50(93) 


0.8300(413) 


[6,12] 


0.95(92) 


0.01296(49) 



TABLE XXI: Decay constants on 12 3 x 24 lattice. Here for the renormalization factor we employ k c determined from a 
simultaneous fit to m 2 PS and m^alh in Table IVI 





fps 




fv 




0.14585 


0.1785(14) 


[6,12] 


0.3118(33) 


[6,12] 


0.14660 


0.15784(87) 


[6,12] 


0.2874(57) 


[6,12] 


0.14705 


0.1413(14) 


[6,12] 


0.2496(97) 


[6,12] 


0.14720 


0.1412(41) 


[6,12] 


0.2422(239) 


[6,12] 



TABLE XXII: Baryon masses on 12 3 x 24 lattice. 





niN 




X 2 /dof 


mA 


\imiri i trnax\ 


X 2 /dof 


0.14585 


1.5357(69) 


[5,12] 


0.65(76) 


1.7722(97) 


[5,12] 


0.74(83) 


0.14660 


1.3619(92) 


[5,12] 


0.85(66) 


1.6061(183) 


[5,12] 


1.45(97) 


0.14705 


1.2054(165) 


[5,12] 


0.69(96) 


1.5110(268) 


[5,12] 


1.28(81) 


0.14720 


1.1791(417) 


[5,12] 


0.99(62) 


1.5300(1020) 


[5,12] 


0.62(1.23) 



TABLE XXIII: Plaquette and rectangular loops on 12 3 x 24 lattice. 





(w lxl ) 


(w lx2 ) 


0.14585 


0.504529(56) 


0.249916(70) 


0.14660 


0.508445(69) 


0.254866(88) 


0.14705 


0.511202(68) 


0.258350(86) 


0.14720 


0.512632(144) 


0.260157(186) 
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TABLE XXIV: Meson and bare AWI quark masses on 16 3 x 24 lattice. 



ftsea 


nips 




X 2 /dof 


my 




X 2 /dof 


m AWI 

' ' L quark 


0.14585 


0.6333(19) 


[6,12] 


0.72(52) 


1.0488(43) 


[6,12] 


0.82(76) 


0.06378(47) 


0.14660 


0.4781(16) 


[6,12] 


3.55(2.04) 


0.9403(70) 


[6,12] 


1.41(93) 


0.03642(40) 



TABLE XXV: Decay constants on 16 3 x 24 lattice. Here for the renormalization factor we employ n c determined from a 
simultaneous fit to m\g and m A ^' k in Table Ivl 



ftsea 


fps 




fv 




0.14585 


0.1804(23) 


[6,12] 


0.3151(45) 


[6,12] 


0.14660 


0.1592(16) 


[6,12] 


0.2913(48) 


[6,12] 



TABLE XXVI: Baryon masses on 16 3 x 24 lattice. 



ft sea 


m N 


\t rain •> ^mtix] 


X 2 /dof 


rriA 




X 2 /dof 


0.14585 


1.5567(91) 


[5,12] 


1.97(92) 


1.7804(113) 


[5,12] 


0.64(54) 


0.14660 


1.3257(118) 


[5,12] 


1.58(91) 


1.5899(124) 


[5,12] 


0.96(77) 



TABLE XXVII: Plaquette and rectangular loops on 16 3 x 24 lattice. 



ft-sea 


(W 1 * 1 ) 


(W lx2 ) 


0.14585 


0.504482(75) 


0.249850(90) 


0.14660 


0.508338(61) 


0.254739(76) 
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FIG. 1: Comparison of convergence with various algorithms for inversions of the quark matrix at n aea = 0.14585 (mps/mv = 
0.60) and n va i ence — 0.14850 {mps/mv = 0.27) on 12 3 x 24 lattice. Conventional BiCGStab corresponds to L = 1. We 
define an iteration Ni ter as a dimension of the Krylov subspace to which approximate solutions belong (3S)l . The number of 
matrix- vector products to obtain an approximate solution is 2 x Nu er - 



12 X24, Ksea = 0.14585 



ii 

1 96 



E 91 



J? 90 ■ 



DS-L^ 

L=4 ^-rj-i 



FIG. 2: The number of iterations required with BiCGStab(L) and BiCGStab(DS-L) for inversions of the quark matrix at 
Ksea = Kvaience ~ 0.14585 (mps/mv ~ 0.60) on 12 3 x 24 lattice. Conventional BiCGStab corresponds to L = 1. 
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FIG. 3: Reversibility violation at large sea quark mass of n sea — 0.14585 (mps/mv = 0.60) (left panel) and small sea quark 
mass of n sea = 0.14705 (m P s/m v = 0.40) (right panel) on 12 3 x 24 lattice. 




Nmolec dt Nmolec dt 

FIG. 4: Effect of the molecular dynamics step size dt on the appearance of spikes in dH = H tria i — Ho at K aea — 0.14660 
(mps/mv = 0.50) on 12' x 24 lattice (left panel). The right panel is an enlargement around the spikes in the case of dt = 0.0025. 
jjD _1 (D t ) _1 0|| with the Wilson-clover operator D and the pseudofermion field <j> as well as the corresponding contribution 
with the smallest eigenvalue Ao and its overlap Co = (xo, 4>) for 75D are also plotted, where Xi is an eigenfunction of 75-D such 
that 4> — CiXi. 
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FIG. 5: Effective masses of pseudoscalar (left panel) and vector meson (right panel) at n S ea = 0.14585 (mps/mv = 0.60) on 
12 3 x 24 lattice. 
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FIG. 6: Effective masses of pseudoscalar (left panel) and vector meson (right panel) at n sea = 0.14705 (mps/mv = 0.40) on 
12 3 x 24 lattice. 
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FIG. 7: Effective potential energies V e //(r = N s /4,t) at K sea = 0.14705 (m P s/m v = 0.40) on 12 3 x 24 lattice. 




FIG. 8: Static quark potentials at K sea = 0.14585, 0.14660, 0.14705 and 0.14720 correspond to m PS /m v = 0.60, 0.50, 0.40 and 
0.35 on 12 3 x 24 lattice. 
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FIG. 10: Sea quark mass dependence of the cumulative autocorrelation time of plaquette on 12 3 x 24 lattice. Open symbols 
are the results obtained in our previous study |3- 




FIG. 11: Bin size dependence of jack-knife error of pseudoscalar meson mass (left panel), and plaquette and rectangular loop 
(right panel) at K sea = 0.14705 {m PS /m v = 0.40) on 12 3 x 24 lattice. 
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FIG. 12: Volume dependence ol pseudoscalar (left panel) and vector meson masses (right panel) at K sea = 0.14585 (mps/mv 
0.60) and K aea = 0.14660 (m PS /m v = 0.50). 
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FIG. 13: Volume dependence ol AWI quark masses at K sea = 0.14585 (mps/mv = 0.60) and K, sea = 0.14660 (mps/mv = 0.50). 
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FIG. 14: Volume dependence of octet (left panel) and decuplet baryon masses (right panel) at n sea — 0.14585 (mps/mv = 0.60) 
and Ksea = 0.14660 (mps/mv = 0.50). 



2.20 
2.15 
2.10 
2.05 
£ 2.00 
1.95 
1.90 
1.85 
1.80 



Ksea = 0.14585 



2.15 



2.05 - 



Ksea = 0.14660 




( 


( 

> 


> 









10 11 12 13 14 15 16 17 18 

N s 



10 11 12 13 14 15 16 17 18 

N s 



FIG. 15: Volume dependence of Sommer scales at K sea — 0.14585 (mps/mv = 0.60) and n sea — 0.14660 (mps/mv = 0.50). 
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FIG. 16: Chiral extrapolation of pseudoscalar meson mass. Open symbols show the results obtained in the previous calculation 
and filled symbols are our new results. Lines are polynomial fits as explained in the figure. 
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FIG. 17: Dependence of % 2 /dof on the fitting range and order of the fitting polynomial (left panel), and that of the critical 
hopping parameter (right panel) for pseudoscalar meson mass. Symbols are placed at the value of l/n sea corresponding to the 
lowest fitting range, which is changed as shown in the figure, while the highest is fixed to n 3sa = 0.1409 (mps/rnv = 0.80). 
Quadratic, cubic and quartic forms as a function of VWI quark mass are tested. 
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FIG. 19: Dependence of X 2 /dof on the fitting range and order of the fitting polynomial (left panel) , and that of the critical 
hopping parameter (right panel) for AWI quark mass. Symbols are placed at the value of 1/nsea corresponding to the lowest 
fitting range, which is changed as shown in the figure, while the highest is fixed to K sea — 0.1409 (mps/mv = 0.80). Quadratic, 
cubic and quartic forms as a function of VWI quark mass are tested. 




FIG. 20: Chiral extrapolation of vector meson mass in terms of pseudoscalar meson mass. Open symbols show the results 
obtained in the previous calculation 3]. 
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FIG. 21: Ratio of vector meson correlators with momentum 2n/L and the polarization parallel and perpendicular to it. 
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FIG. 22: Chiral extrapolation of pseudoscalar (left panel) and vector (right panel) meson decay constants. Open symbols show 
the results obtained in the previous calculation Q. 
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FIG. 23: Chiral extrapolation of octet (left panel) and decuplet (right panel) baryon masses. Open symbols are the results in 
our previous study Q. 
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FIG. 24: Chiral extrapolation of ro. Open symbols are the results in our previous study Q- 
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FIG. 25: Comparison of degenerate up and down quark masses obtained by chiral extrapolations with polynomials. Open 
symbols show the results obtained in the previous calculation 3] and filled symbols are our new results. Lines are combined 
linear continuum extrapolations in the previous calculation. 




FIG. 26: Test of simultaneous continuum ChPT fit to pseudoscalar meson mass and decay constant. In this plot, quark mass 
defined through the axial vector Ward identity is used. The right panel shows the ratio m% s /2rn^^' k to focus on the chiral 
logarithm behavior. Open symbols are the results obtained in our previous study 
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FIG. 27: Test of simultaneous continuum ChPT fit with the quark mass defined through the vector Ward identity. Open 
symbols are the results obtained in our previous study 




FIG. 28: Chiral extrapolation of vector meson mass with a polynomial in Eq. 134H and a function motivated by ChPT in 
Eq. 1 12t . Open symbols are the results obtained in our previous study Q. 




FIG. 29: Chiral extrapolation of octet (left panel) and decuplet (right panel) baryon masses with polynomials in Eq. 1361 1 and 
functions motivated by ChPT in Eq. H43II . Open symbols are the results obtained in our previous study 0. 
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FIG. 30: Test of the WChPT fit to pseudoscalar meson mass, AWI quark mass and decay constant. The right panel shows the 
ratio m 2 PS /2m^^' k to focus on the chiral logarithm behavior. Open symbols are the results obtained in our previous study Q. 




FIG. 31: Test of the resummed WChPT fit to pseudoscalar meson mass and AWI quark mass. The right panel shows the ratio 
mp S /2m^' t to focus on the chiral logarithm behavior. Open symbols are the results obtained in our previous study 0]. 




FIG. 32: Ratio of the next-to-leading order term to the leading one for m PS with the resummed WChPT formulae (left panel) 
and with WChPT formulae without resummation (right panel) as a function of mYuark- 
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FIG. 33: Comparisons of the polynomial and the resummed WChPT fits to pseudoscalar meson mass and AWI quark mass 
determined at mps/mv = 0.80-0.35. Circles show the lattice data and the square is the extrapolated result at the physical 
point. Open symbols are the results obtained in our previous study Q- 




FIG. 34: Comparison of quadratic and resummed WChPT fits to pseudoscalar meson masses and AWI quark masses determined 
from the previous data of rnps/mv = 0.80-0.55 (open symbols) with the new small sea quark mass data (filled symbols). 
The right panel is an enlargement around the chiral limit. 




FIG. 35: Continuum extrapolations of degenerate up and down quark mass obtained by chiral extrapolations with polynomials 
Q (open symbols) and the resummed WChPT formulae (filled symbols). The star at a ~ lGeV -1 ((3 — 1.8) represents the 
results obtained by the resummed WChPT formulae with data at mps/rriv ~ 0.80-0.35. The others are the results with 
mps/mv = 0.80-0.55. The dashed lines are the combined linear fit to the quadratic chiral fit results and the dashed-dot lines 
are the ones to the resummed WChPT fit results, both with mps/mv = 0.80-0.55. The solid lines are the combined linear fits 
to the resummed WChPT chiral fit results with our whole data of mps/mv = 0.80-0.35 at ft = 1.8 and mps/mv = 0.80-0.55 
at (3 = 1.95 and 2.1. 



